perm filename LOCUS.SAI[1,BGB] blob sn#001248 filedate 1972-10-22 generic text, type T, neo UTF8
00100	ENTRY LOCUS;
00200	BEGIN	"LOCUS"
00300		REQUIRE "ABBREV[SYS,BGB]" SOURCE_FILE;
00400		REQUIRE "COMMON[IA,BGB]" SOURCE_FILE;
00500		REQUIRE "TRIGER[SYS,BGB]" SOURCE_FILE;
00600	EXTERNAL INTEGER PROCEDURE LS1V3P (REAL ARRAY V,ARGS);
00700		REAL ARRAY ITEMVAR CAM,F;
00800		INTEGER M,N,FLG,CLNCNT;
00900		SET FSET,LTSET;
01000		REAL ARRAY ∂CAM[1:4,1:3];
01100		REAL ARRAY V1,V2,V3[1:3];
01200		REAL ARRAY V[1:6,1:3];
01300		REAL ARRAY ARGS[1:4,1:3];
01400		REAL ARRAY ∂LT[0:3];
01500		DEFINE X(W)="∂(W)[1]",Y(W)="∂(W)[2]",Z(W)="∂(W)[3]";
     

00100	PROCEDURE GET;
00200	BEGIN	"GET"
00300		SET CAMSET;
00400	α GET OR CREATE A CAMERA FOR THE TVFILE;
00500		CAMSET	←	LOCOR⊗TVFILE;
00600		IF LENGTH(CAMSET)=0 THEN
00700	BEGIN
00800		OUTSTR("NEW CAMERA MADE"&↓);
00900		CAM	←	NEW(∂CAM);
01000		MAKE LOCOR⊗TVFILE≡CAM;
01100	END	ELSE
01200		CAM	←	LOP(CAMSET);
01300	α REPORT RASTER SIZE AND FOCAL;
01400		SETFORMAT(0,4);
01500		OUTSTR("	FOCAL	"&CVG(FOCAL*304.8)&" mm"&↓);
01600		OUTSTR("	RASTER SIZE	");
01700		OUTSTR(CVG(2*PDY*304.8)&" mm  by  ");
01800		OUTSTR(CVG(2*PDX*304.8)&" mm"&↓);
01900		SETFORMAT(0,7);
02000	α GET ALL THE FEATURES IN THIS PICTURE;
02100		∀ F|FεFEASET ∧ TVFILE⊗F≡ANY DO PUT F IN FSET;
02200		N	←	LENGTH(FSET);
02300		M	←	N*(N-1)*(N-2)/6;
02400		OUTSTR(↓);
02500		OUTSTR(9&CVS(N)&" LANDMARK POINTS."&↓);
02600		OUTSTR(9&CVS(M)&" LANDMARK TRIANGLES."&↓);
02700	END	"GET";
     

00100	α CREATE LANDMARK TRIANGLES;
00200	PROCEDURE SETUP;
00300	BEGIN	"SETUP"
00400		REAL ARRAY ITEMVAR F1,F2,F3,WN1,WN2,WN3,LT;
00500		∀ F1,F2,F3|F1εFSET∧F2εFSET∧F3εFSET DO
00600	BEGIN	"LOOP"
00700		LABEL EOL;
00800	α FETCH FEATURE WINDOWS;
00900		IF #(F1)<#(F2) ∧ #(F2)<#(F3) THEN ELSE GO EOL;
01000		WN1	←	COP(TVFILE⊗F1);
01100		WN2	←	COP(TVFILE⊗F2);
01200		WN3	←	COP(TVFILE⊗F3);
01300	α INSURE CLOCKWISE ORDER;
01400	BEGIN	"CWORDER"
01500		REAL A,B,C,Q,R;
01600		A	←	Y(WN1) - Y(WN2);
01700		B	←	X(WN2) - X(WN1);
01800		C	←	X(WN1)*Y(WN2) - X(WN2)*Y(WN1);
01900		Q	←	A*X(WN3) + B*Y(WN3) + C;
02000		IF Q>0 THEN BEGIN F2↔F3;WN2↔WN3;END;
02100	α COLINEAR TEST  -  IN IMAGE PIXELS;
02200		R	←	ABS(Q / SQRT(A↑2 + B↑2));
02300		IF R<50 THEN BEGIN CLNCNT←CLNCNT+1;GO EOL;END;
02400	END	"CWORDER";
02500	
02600	α THE COEFFICIENTS OF A PLANE;
02700	BEGIN	"COEF"
02800		REAL AX,AY,AZ,BX,BY,BZ,K;
02900		AX  ←  X(F1)-X(F2);		BX  ←  X(F2)-X(F3);
03000		AY  ←  Y(F1)-Y(F2);		BY  ←  Y(F2)-Y(F3);
03100		AZ  ←  Z(F1)-Z(F2);		BZ  ←  Z(F2)-Z(F3);
03200		∂LT[1] ← (AY*BZ - AZ*BY);
03300		∂LT[2] ← (AZ*BX - AX*BZ);
03400		∂LT[3] ← (AX*BY - AY*BX);
03500		K	←	∂LT[1]*X(F1) + ∂LT[2]*Y(F1) + ∂LT[3]*Z(F1);
03600		∂LT[1]←∂LT[1]/K;
03700		∂LT[2]←∂LT[2]/K;
03800		∂LT[3]←∂LT[3]/K;
03900	END	"COEF";
04000	
04100	α CREATE A LANDMARK TRIANGLE;
04200		LT	←	NEW(∂LT);
04300		PUT LT IN LTSET;
04400	α PUT ORDERED FEATURE ITEM NUMBERS INTO LT'S DATUM;
04500	BEGIN	"PUT"
04600		INTEGER ARRAY ITEMVAR ILT;
04700		ILT	←	LT;
04800		∂(ILT)[0]←	(#(F1)ROT -12) LOR (#(F2)LSH 12) LOR #(F3);
04900	END	"PUT";
05000	EOL:
05100	END	"LOOP";
05150		OUTSTR(9&CVS(LENGTH(LTSET))&" LANDMARK TRIANGLES THRU SETUP."&↓&↓);
05200	END	"SETUP";
     

00100	PROCEDURE GETCOSINES (INTEGER ARRAY ITEMVAR ILT);
00200	BEGIN	"GETCOSINES"
00300		DEFINE THRICE="FOR I←1 STEP 1 UNTIL 3 DO";
00400		REAL ARRAY ITEMVAR F1,F2,F3,WN1,WN2,WN3;
00500		INTEGER NNN,I;	REAL R;
00600	α GET THE FEATURES AND THEIR WINDOWS FOR THIS TRIANGLE;
00700		NNN	←	∂(ILT)[0];
00800		F1	←	CVI((NNN ROT  12) LAND '7777);
00900		F2	←	CVI((NNN ROT -12) LAND '7777);
01000		F3	←	CVI(NNN LAND '7777);
01100		WN1	←	COP(TVFILE⊗F1);
01200		WN2	←	COP(TVFILE⊗F2);
01300		WN3	←	COP(TVFILE⊗F3);
01400	α STASH REAL WORLD 3D FEATURE LOCUS;
01500		FOR NNN←1 STEP 1 UNTIL 3 DO THRICE
01600		ARGS[NNN,I]←(CASE(NNN-1)OF(∂(F1)[I],∂(F2)[I],∂(F3)[I]));
01700	α COMPUTE THE COSINES BETWEEN RAYS TO FEATURES;
01800		V1[1]	←	PDX*X(WN1)/LDX;
01900		V1[2]	←	PDY*Y(WN1)/LDY;
02000		V1[3]	←	-FOCAL;
02100		R	←	SQRT(V1[1]↑2 + V1[2]↑2 + V1[3]↑2);
02200		THRICE	V1[I] ← V1[I] / R;
02300	
02400		V2[1]	←	PDX*X(WN2)/LDX;
02500		V2[2]	←	PDY*Y(WN2)/LDY;
02600		V2[3]	←	-FOCAL;
02700		R	←	SQRT(V2[1]↑2 + V2[2]↑2 + V2[3]↑2);
02800		THRICE	V2[I] ← V2[I] / R;
02900	
03000		V3[1]	←	PDX*X(WN3)/LDX;
03100		V3[2]	←	PDY*Y(WN3)/LDY;
03200		V3[3]	←	-FOCAL;
03300		R	←	SQRT(V3[1]↑2 + V3[2]↑2 + V3[3]↑2);
03400		THRICE	V3[I] ← V3[I] / R;
03500	
03600		DEFINE DOT(X,Y)="X[1]*Y[1]+X[2]*Y[2]+X[3]*Y[3]";
03610		ARGS[4,1]←	DOT(V2,V3);
03620		ARGS[4,2]←	DOT(V1,V3);
03630		ARGS[4,3]←	DOT(V1,V2);
03640	
04000	END	"GETCOSINES";
     

00100	α RUN ALL THE LANDMARK TRIANGLES THRU THE LOCUS SOLVER;
00200		PROCEDURE RUN;
00300	BEGIN	"RUN"
00400		REAL XX,YY,ZZ;
00500		INTEGER I,J,NOCNT,WINCNT,DOUBLES,NROOTS,WINNER,HLFCNT;
00550		REAL ARRAY ITEMVAR LT,LT2;
00600	α INITIALIZATION;
00700		XX←YY←ZZ←0;
00800		HLFCNT←NOCNT←WINCNT←DOUBLES←0;
00900		∀ LT|LTεLTSET DO
01000	BEGIN	"SOLVE"
01100		GETCOSINES(LT);
01200		NROOTS ← LS1V3P(V,ARGS);
01300		IF NROOTS≤0 THEN NOCNT←NOCNT+1;
01400		FOR I←1 STEP 1 UNTIL NROOTS DO
01450	BEGIN
01500		FOR J←1 STEP 1 UNTIL 3 DO 
01600		OUTSTR(9&CVG(V[I,J]));OUTSTR(↓);END;
     

01400	α A WINNER IS ONE SOLUTION IN ALL THE PROPER HALFSPACES;
01600		FOR J←1 STEP 1 UNTIL NROOTS DO 
01700	BEGIN	"CHECKER"
01800		DEFINE X="V[J,1]",Y="V[J,2]",Z="V[J,3]";
02000		WINNER ← TRUE;
02100		∀ LT2|LT2εLTSET DO
02200		IF ∂(LT2)[1]*X + ∂(LT2)[2]*Y + ∂(LT2)[3]*Z < 1 THEN
02300		BEGIN WINNER←FALSE;OUTSTR("X"&↓);HLFCNT←HLFCNT+1;DONE;END
02350		ELSE OUTCHR(".");
02400		IF WINNER THEN BEGIN WINNER ← J;OUTCHR("!"&13&10);DONE;END;
02500	END	"CHECKER";
02600		IF WINNER THEN
02700	BEGIN	"ACCUMULATE"
02800		FOR I←1 STEP 1 UNTIL 3 DO OUTSTR(9&CVG(V[WINNER,I]));OUTSTR(↓);
02900		XX←XX+V[WINNER,1];
03000		YY←YY+V[WINNER,2];
03100		ZZ←ZZ+V[WINNER,3];
03200		WINCNT←WINCNT+1;
03300	END	"ACCUMULATE";
03333	END	"SOLVE";
03400		OUTSTR(9&CVS(WINCNT)&" WINNERS"&↓);
03500		OUTSTR(9&CVS(NOCNT+CLNCNT+DOUBLES)&" LOSERS"&↓&↓);
03600		OUTSTR(9&CVS(CLNCNT)&" COLINEARS"&↓);
03700		OUTSTR(9&CVS(NOCNT)&" NO ROOTS"&↓);
03800		OUTSTR(9&CVS(DOUBLES)&" DOUBLE ROOTS"&↓);
03900		OUTSTR(↓);
04000		OUTSTR(9&CVG(XX/WINCNT));
04100		OUTSTR(9&CVG(YY/WINCNT));
04200		OUTSTR(9&CVG(ZZ/WINCNT));
04400	END	"RUN";
     

00100	INTERNAL PROCEDURE LOCUS;
00200	BEGIN	"LOCUS"
00300		GET;
00400		SETUP;
00500		RUN;
00600		OUTSTR(↓&"*");
00700	END	"LOCUS";
00800	END	"LOCUS";